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Abstract 

We present a new Fortran library to evaluate all harmonic polylogarithms up to 
weight four numerically for any complex argument. The algorithm is based on a 
reduction of harmonic polylogarithms up to weight four to a minimal set of basis 
functions that are computed numerically using series expansions allowing for fast 
and reliable numerical results. 
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1 Introduction 



Feynman integrals in perturbative quantum field theory are generically ex- 
pressed in terms of the classical polylogarithm functions Li n (z) and the Nielsen 
polylogarithms S n ^(z) [I]. In the late nineties, it was realized that these classes 
of functions are too restricted when going beyond one-loop level in the pertur- 
bative expansion, where new functions appear that can no longer be expressed 
in terms of the classical polylogarithm functions. While a completely generic 
generalization of polylogarithms has been studied in the mathematical liter- 
ature (going under the name of multiple polylogarithms [2113] ). it is mostly 
only a specific subset of multiple polylogarithms, the so-called harmonic poly- 
logarithms [1] and their two-dimensional and cyclotomic generalizations 06], 
that make their appearance in the theoretical predictions of physical quantities 
beyond leading-order. 

In this paper we concentrate exclusively on harmonic polylogarithms (HPL's) 
up to weight four, which since their introduction have found many applica- 
tions in computations up to two-loop order in the perturbative expansion, 

der to confront the theoretical next-to-next-to-leading order (NNLO) pre- 
dictions to experiment, it is mandatory to be able to evaluate HPL's nu- 
merically in a fast and accurate way. The requirements to such a numeri- 
cal code are twofold: first, the evaluation should be fast, because the use of 
NNLO matrix elements in Monte Carlo integration codes may require thou- 
sands, if not millions, of function calls. Second, it is desirable to be able to 
compute HPL's for arbitrary complex arguments, which appear for example 
when the complex mass scheme is employed or for certain kinematic config- 
urations in loop calculations involving massive particles [T^PU] . In the last 
decade, various codes have been developed to evaluate HPL's numerically. 
While the code hplog [30J, written in Fortran, is restricted to the evalua- 
tion of HPL's up to weight four and for real values of the arguments, the 
code HPL (Mathematica)[31,32| and the implementation of the harmonic 
polylogarithms into the GiNaC framework (C++) [33] are generic and allow to 
evaluate in principle any harmonic polylogarithm with arbitrary precision for 
any complex argument. 

The focus of this paper is the Chaplin (Complex HArmonic PolyLogarithms 
In fortraN) library, a new Fortran code that allows to evaluate numerically 
all harmonic polylogarithms up to weight four for arbitrary complex argu- 
ments. While Chaplin is similar in spirit to the aforementioned codes, its 
main advantages lie in speed, through the use of Fortran as a programming 
language, and in its capability to compute HPL's numerically for any point 
in the complex plane. Chaplin reduces each of the 120 HPL's up to weight 
four to a set of 32 basis functions |36j, which are entirely expressed through 
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only two new functions of weight four besides the classical polylogarithms. 
The basis functions are then mapped to the interior of the unit circle, where 
they are computed numerically using suitably chosen series expansions that 
allow to obtain a fast numerical convergence. 

The paper is organized as follows: In Section [2] we give a short review of 
harmonic polylogarithms and of their main algebraic and analytic properties. 
In Section [3] we review the reduction of all HPL's up to weight four to the 
set of basis functions introduced in Ref. [36]. The series expansions used by 
Chaplin to compute the basis functions numerically are derived in Section HJ 
while the Chaplin library itself, together with comparisons to hplog, HPL 
and GiNaC, is presented in Section 



2 Short review of harmonic polylogarithms 



In this section we give a short review of harmonic polylogarithms (HPL's), as 
they are at the heart of the Chaplin library. HPL's are defined recursively 
via the iterated integrals [I] 

H(a 1 ,...,a n ;z) = [ dt f(ar,t) H(a 2 , . . . ,a n ;t) , (2.1) 
J o 

where aj G {—1, 0, 1} and 

/(-M) = TT7' /(M = 7> f0-' z ) = v~ ' ( 2 - 2 ) 

and where we defined H(; z) = 1. If all the aj are simultaneously zero, the 
integral (12 .1}) is divergent, and so in this case we use the definition, 

H(6 n ,z) = ±-log n z, (2.3) 
n\ 

where we used the obvious vector notation a n = (a,..., a). Note that the 

n times 

number n of indices is usually referred to as the weight of the HPL. Harmonic 
polylogarithms are a generalization of the classical polylogarithm functions, 
defined recursively by, 

r z dt 

U n (z)= — Li n _i(£) and Lii(z) = - log(l - z) . (2.4) 
Jo t 



Iterated integrals are well-known to form a shuffle algebra, and so in particular 
we can express a product of two harmonic polylogarithms of weight n\ and n 2 
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clS db linear combination of HPL's of weight n\ + ri2, 

H(a 1 ;z)H{a 2 ;z) = J2 H &z), (2-5) 

a=ai III 3,2 

where a\ HIS2 denotes the shuffle of the two weight vectors a\ and a 2 , i-e., all 
possible concatenations of ai and a 2 in which relative orderings of cti and a 2 
are preserved. 

Up to weight three, HPL's are known to be expressible through ordinary log- 
arithms and the classical polylogarithms Li n only. Starting from weight four, 
not all HPL's can be expressed in terms of classical polylogarithms, and gen- 
uine new functions appear. For special classes of HPL's, however, it is possible 
to find closed expressions in terms of other functions. An example of this was 
already given in Eq. ( I2.3p . Moreover, we have, 

H(s n ;z) = ^-^\og n (l-sz) and H(0 n ^, 1; z) = U n {z) , (2.6) 

i.e., harmonic polylogarithms contain the classical polylogarithms as special 
cases. 

Apart from the shuffle relation (12.51) . HPL's satisfy various intricate func- 
tional equations, relating HPL's with different arguments among each other. 
As an example, the functional equation relating HPL's with opposite argu- 
ments reads, 

H(a;-z) = {-If H{-a;z) , (2.7) 

where p denotes the number elements in a equal to ±1. Furthermore, it is 
always possible to express harmonic polylogarithms of the form H(a; 1/z) as 
a linear combination of HPL's of the form H(a; z). This allows in particular 
to analytically continue the harmonic polylogarithms outside the unit disc, a 
property that will be used later in the numerical implementation of the HPL's 
into the Chaplin library. Additional relations among HPL's with related 
arguments have been presented in Ref. [lf3T] . 

Let us conclude this section by discussing some special values of the argument 
z for which the HPL's can be expressed in terms of known transcendental 
numbers. First, it is easy to see that, unless a = n , all HPL's vanish for 
2 = 0. Second, if the argument z of a harmonic polylogarithm is ±1, then it 
can be expressed in terms of so-called colored multiple zeta valued] (CMZV's), 

#(0 1 _ UI 0, Sl ,...,0 i _^ I 0, Sfe ;l) 

m\—\ m^—l (2.8) 

= (-l) m+P C(mi, ...,m k ; s h s 2 /si, . . . , s fc /s fc _i) , 
1 Also known as alternating multiple zeta values or Euler-Zagier sums. 
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with S{ = ±1 and m = mi + . . . + ntf., and p denotes the number of elements 
in {s^ equal to +1. The CMZV's are defined by nested sums, 

cr ni cr "2 . . . a n k 

C(m l ,...,m k ;a 1 ,...,a k ) = ]T — - 2 m2 , (2.9) 

0<m<n 2 <-<n fc n l U 2 

with (jj = ±1. CMZV's are convergent if and only if (mi,cri) ^ (1>1)- It 
follows then immediately that HPL's of the form H(±l, a; ±1) are in general 
divergent!!]. 



3 Reduction to basis functions 

From the previous section it is clear that many HPL's are not independent 
functions, but they are related among themselves by various intricate relations. 
In order to achieve an efficient numerical implementation, it is desirable to 
have as few independent functions as possible, i.e., we would like to resolve 
all the identities in order to arrive at a minimal set of functions, which are 'as 
simple as possible' and through which all other HPL's can be expressed. 

All the functional equations among harmonic polylogarithms (or, more gener- 
ically, among multiple polylogarithms) can be resolved through the so-called 
symbol calculus. At the heart of the symbol calculus is the so-called symbol 
map [34J, a linear map that associates to an HPL of weight n a tensor of 
rank n. As an example, the tensor associated to the classical polylogarithm 
Li n (z) = if(O n _i, 1; z) reads, 

S(U n (z)) = -(1 - z) <g> z ® . . . g z . (3.1) 

(n— 1) times 

Furthermore, the symbol maps products that appear inside the tensor product 
to a sum of tensors, 

...®(X-Y)®... = ...®X®... + ...®Y®.... (3.2) 

It is conjectured that all the functional identities among (multiple) polyloga- 
rithms are mapped under the symbol map S to algebraic relations among the 
tensors. Hence, the symbol calculus provides an effective way to resolve all the 
functional equations among (a certain class of) multiple polylogarithms. 

In Ref. [36], the symbol map was used to obtain a set of basis functions through 
which all HPL's up to weight four can be expressed. The basis functions ob- 
tained in Ref. [36] read, 

2 In some cases the divergence can be tamed, e.g., lim z _;.o H(l, 0; z) = — £2- 
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for weight one, 

B { l\z) = logz, B?\z)=\og(l-z), B?\z) = \og(l + z), (3.3) 
for weight two, 

JSlp{z) = U*(z), B<£\z) = U 2 {-z), B < i\z) = U 2 ( 1 —^] , (3.4) 



for weight three, 

B^\z) = U 3 (z), Bf\z) = his(-z), Bf\z) = Li 3 (l - z) 

4 4) W = Li 3 ( T ^), B?(z) = U 3 {^), Bf\z) = U, n ~ Z 

B«\z) =uJ\^-\ , Bf\z) = u/ 2Z 



(3.5) 



1 + 2/ ' ° ' 7 

for weight four, 

B<i\z) = Li 4 (z), B?\z) = LU{-z), 

Bf\z) = LU(1 - z), £ 4 4) (^) = Li 4 (^J-) , 

4 5) (*) = Li 4 ( 7 ^ T ) , 4 6) (*) = Li 4 (^) , 

M 7) W = L l4 (^), B f) (z) = Ii4 (Lzf), 

bW(z) = uJ^], B^(z) = U 1 1 



z 2 



Bf\z) = Li 4 (l - , 2 ) , B^\z) = Li 4 



B^\z) = U 4 



Az 



All harmonic polylogarithms up to weight three can be expressed through the 
basis functions in Eq. (13.31 - 13. 6j) . Starting from weight four, we need to extend 
the set of functions by adjoining three new elements to the basis, 

4 16) W = Li 2 , 2 (-M), B[ 17 \z) = U 2 J 1 - 21 



B { A 18) (z) = Li, 



j 4 yt) — uj.2,2 77) 



1 2z 



2'z+lJ' (3J) 



.2' z-l, 

where Li 2 2 denotes a two-variable multiple polylogarithm that cannot be ex- 
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pressed through classical polylogarithms only, 



oo ni-1 n\ z n 2 

Li 2 , 2 (*i, *a) = ^3^' ( 3 - 8 ) 

ni=ln 2 =l 1 ' l 2 

For practical purposes we find it more convenient to use a different set of 
multiple polylogarithms as basis functions than the one used in Ref. |36j . 
More specifically, we find it more convenient to perform a change of basis and 
replace the functions B^\z), for i = 16, 17, 18, by the functions B^\z), which 
are directly expressed as HPL's, 

B^\z) =H(0,l,0,-l;z) = B { 4 16 \-z), 

BT\z) =#(0,1,1,-1;*), (3.9) 

B^\z) =#(0,1,1,-1;-*). 

The set of the 32 functions B^\z) defines a basis through which all HPL's 
up to weight four can be expressed. As a consequence, any numerical code to 
evaluate this set of basis functions will automatically be able to evaluate all 
120 HPL's up to weight four. Furthermore, as the basis functions only involve 
the two genuine multiple polylogarithms #(0, 1, 0, — 1; z) and #(0, 1, 1, — 1; z) 
besides the ordinary logarithms and the classical polylogarithms Li 2 (*), Li 3 (z) 
and Li 4 (z), it is enough to have numerical routines for these latter functions. 
In this way we can reduce the problem of evaluating the 120 HPL's up to 
weight four to only a handful of non-trivial numerical routines. In the CHAP- 
LIN library, these routines consist of series expansions for the aforementioned 
functions that will be described in the next section. 

Let us conclude this section by reviewing some of the properties of the basis 
functions B^\z) derived in Ref. [36]. First, it is easy to check that all the 
basis functions are real for z G [0, 1]. However, we stress that the expressions 
in Eq. (13.31 - 13. 9p are strictly valid only for z G [0,1]. While most of the 
expressions are valid everywhere throughout the unit disc, the analytic form 
of B± (z) valid on the whole interior of the disc reads (36], 

Bf\z) 

' Li 4 (l - z 2 ) , if Re(z) > or (Re(z) = and Im(z) > 0) , 

Li 4 (l — z 2 ) — y a(z) log 3 (l — z 2 ) , otherwise , 

(3.10) 

where a(z) = sign(Im(*)). Second, since it is our goal to build a numerical 
code to evaluate HPL's for arbitrary complex arguments, we need to analyti- 
cally continue the basis functions outside the unit disc. In Ref. [36] inversion 
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relations of the form 



B^\z) = Cijki E>k (~ ) + products of lower weight. (3.11) 



3 

k,l 



were derived that can be used for this purpose. Finally, we note that there is 
a subtlety in the basis function B^ (z) when going from the interior to the 
exterior of the unit disc because B^ (z) has a branch cut along the unit circle 
in the complex z-plane. In Ref. f36] it was shown that if we want B^ 5 \z) to be 
continuous and real for z G [0, 1], we need to choose the following prescription 
for \z\ = 1, 

) (z) = Li 4 + iv(*)e) ■ (3.12) 



4 Numerical evaluation of the basis functions 



4-1 Notations and conventions 



In the previous section we introduced a set of basis functions through which 
every HPL up to weight four can be expressed. The basis has the property 
that it only involves two types of new functions, besides the ordinary loga- 
rithm and the classical polylogarithms. These two new functions can be cho- 
sen to correspond to the two harmonic polylogarithms H(0, 1, 0, — 1; z) and 
H(0, 1, 1, — 1; z). In this section we present series expansions of these func- 
tions which are used inside Chaplin to evaluate the basis functions. 

Let us start by introducing some notations and conventions. As we will deal 
with series expansions, we define some operations on the coefficients of the 
series, i.e., on sequences of complex numbers. For two sequences of complex 
numbers a n and b n , n G N, we define their convolution product as the sequence 
(a * b) n defined by 

(a*b) n = J2 \i ) a kK-k- (4.1) 
k% W 

It is easy to check by manipulating the sum that this operation is associative, 
commutative and has the sequence e n = 5o )n , where 5{j is the Kronecker 
symbol, neutral element, 

a * (b * c) = (a * b) * c , a * b = b * a , a*e = e*a = a. (4-2) 

The fact that e n is a neutral element is obvious, and the commutativity follows 
immediately from changing the summation variable from k to n — k. Associa- 
tivity is less obvious, and is proved in Appendix[A] Furthermore, this operation 
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is compatible with the usual termwise addition and scalar multiplication of 

sequences, 

a*(b + c) = a*b + a*c, (k • a) * b — k ■ (a * b) , (4.3) 

where a, b and c are sequences of complex numbers and k is a constant com- 
plex number. The convolution product allows us to write the coefficients that 
appear in the product of two power series as the convolution of the coefficients 
of the individual factors, e.g., 



oo 



E^ m ) (eV] = E^#^- (4-4) 



Finally, for later convenience, we define for a given sequence a n of complex 
numbers the three new sequences a n , a n and s(a) n by 

, a n -\ , if n > 1 
-l) n a n , s(a) n ={ ~ . (4.5) 



n + 1 I , otherwise 

The bar-operation allows us to define the coefficients of the series expansion 
of f(—z) in terms of the coefficients of the series expansion of f(z). More 
precisely, the two series expansions are related by 



oo 



/(*) = E fn z n and /(-*) = J2fnZ n - (4.6) 

n=0 n=0 

We also define the composition of the 'o '-operation and the convolution and 
shift operations, 

( I U\o ( a * b )n , o, x s(a) n . . 

(a © b) n = (a * 6) n = — and s{a) n = — — . (4.7) 

n+l n+l 

Note that the ©-operation is commutative, but not associative. For later con- 
venience we introduce the following convention, 

a®b®c = a®(b®c). (4.8) 

The operations on sequences of complex numbers we just defined allow us to 
write the coefficients appearing in the series expansion of the basis functions in 
a compact closed form. The sequences of complex numbers that appear inside 
these closed-form expressions are well-known sequences of (rational) numbers 
which we recall in the following. 

(1) (Shifted) C values: 

C fc) = { (4.9) 
if fe _i, if k - n = 1 , 
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where ( m = C( m ) denotes the Riemann ( function and H m the m-th 
harmonic number, 

OO 1 TO 1 

CW = Ezi and tf m =£-. (4-10) 

n=l '* n=l 

Note that is rational for n > k and transcendental of weight k — n 
otherwise. 

(2) Bernoulli numbers: The Bernoulli numbers B n are defined through the 
generating series 

~ oo n 

T £ T = E S » £ T- 

Note that B 2n +i = 0, Vn G N. The Bernoulli numbers are related to ( 
values via the formula, 

Co = #i = ~\ and C-n = for ra > 1 . (4.12) 

(3) Genocchi numbers: The Genocchi numbers are defined through the 
generating series 

o- oo fi 

They are related to the Bernoulli numbers via 

G n = 2(1 -T)B n . (4.14) 

Note that G = G 2n +i = 0. 

(4) Polylogarithms in half-integer values: The series expansion of the 
basis functions also involve the following sequences of numbers 

4 = (-l)» Li_ n Q . (4.15) 

The sequence £ n admits the generating function, 

1 -EEL 7 n 

=r=S*-sf- (4 - 16) 



Ze " n=0 



^.2 Numerical evaluation of classical polylogarithms 



In this section we discuss series expansions of classical polylogarithms that 
can be used to obtain reliable numerical results for these functions (at least 
in some regions of the complex plane). In particular, truncated versions of 
these series are used by Chaplin to evaluate the classical polylogarithms. 
As the goal of CHAPLIN is to provide an efficient way to evaluate harmonic 
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polylogarithms for arbitrary complex arguments, we divide the problem into 
two regions: for complex numbers z with \z\ < 1 we use the series expansions 
to evaluate the basis functions numerically, whereas points with \z\ > 1 are 
mapped back into the interior of the unit disc using the inversion formulae for 
the classical polylogarithms. Hence, from here on we will only concentrate on 
complex number z with \z\ < 1. 

Classical polylogarithms can be expanded into a power series around z = 0, 

Lin(*) = £4 (4-17) 
fe=i K 

Even though this series is convergent for \z\ < 1, the convergence is rather 
slow. A faster convergence can be achieved by using the so-called Bernoulli 
substitution [33J, which consists in expanding Li n (z) into a series in log(l — z). 
While this expansion converges much faster for \z\ *C 1 than the Taylor ex- 
pansion (I4.17p . it fails to produce reliable results when z approaches 1. In 
Ref. [35J an alternative expansion of the classical polylogarithms into a series 
in \ogz was derived. In this case the convergence is fast inside an annulus 
around z = 0, but fails to converge for \z\ <C 1. The strategy seems thus clear: 
we can split the interior of the unit disc into two distinct regions, and in each 
region one of the two series expansions converges quickly. Similar expansions 
can also be derived for the two remaining basis functions, H (0, 1, 0, —1; z) and 
H(0, 1, 1, —1; z) (and in principle for every HPL) and are discussed in the rest 
of this section. We start by deriving in detail the expansions of the diloga- 
rithm, because even though these results are well-known, the techniques used 
in the derivation will be the starting point for the higher-weight cases. 



Series expansions of Li 2 (z). Let us start by deriving the expansion of 
Li 2 (z) into a series in log(l — z). Letting x = — log(l — z), this is equivalent 
to finding the Taylor series expansion of the function Li 2 (l — e~ x ). We start 
from the integral representation f!2.4[) and we get, 

Li 2 (l - e~») = / -Lhit) = - - log(l - t) . (4.18) 

Jo t Jo t 

Performing the change of variables t = 1 — e~*' and using Eq. ( 14. lip , we obtain, 

(4.19) 

or equivalently in terms of the original variable z, 

00 r> 

Li 2 (z) = E 77-^ (- Ml - ■ (4.20) 
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The series expansion (14.201) converges rather quickly inside a disc around z = 
of radius R < 1 (the precise value of R used by Chaplin will be given in the 
next section). In the remaining annulus R < \z\ < 1 the dilogarithm admits a 
series expansion in log z 



(4.21) 



oo a(2) 

u 2 (z) = - log z iog(- log z) + -fr logfc z 



Let us sketch the derivation of Eq. (14.211) . Letting x = \ogz, we start from the 
integral representation of the dilogarithm and perform the change of variable 



t = e* 



re" r\f r-x 

U 2 (e x ) =( 2 + - Lh(t) = ( 2 + dt' Li^e* ) . (4.22) 
Ji t Jo 

In order to proceed, we need the Taylor expansion of Li 1 (e x ) = — log(l — e x ). 

Using the integral representation of Lii as well as Eq. (14. lip , we obtain, 



eX dt 



re 

Li ^> = /„ T 



lim 

t e^O 



lim 



log(l-e £ )- 

rx Af °° R 

E 

n=l 



dt' (-If) 



V e-*' - 1 



(4.23) 



The last term in Eq. ( 14 . 2 3 [) is finite, whereas the logarithmic divergences cancel 
between the first two terms, 



lim 



logfl — e £ 



r x dt'~ 






= lim 


L T _ 





log(-e + 0{e 2 )) - log(-x) + log(-e) 



lim [— log(l + 0(e)) — log(— x)\ = — log(— x) . 



(4.24) 



Hence, using Eq. (14. 12j) and the fact that £_ n = for even n, we get, 

C— n 



00 R (— r ) n 
Lii(e^) =-log(-x)-E ' 



n=l 
00 



72! n 



•log(-x) + Cox + J2 



i(n + l)! 



(-*) 



n+l 



l0S( -* )+ £(^V" +1 



(4.25) 



Inserting Eq. (14. 25ft into Eq. ( 14. 22ft and integrating term by term immediately 
reproduces Eq. ( 14 . 2 1|) (after identification x = log z). Note the appearance of 
the nested logarithm, log(— logz) in Eq. (j4.21j) . which seems to be divergent 
for z close to 1. It is however easy to check that 



lim log z loe(— log z) = 



(4.26) 



and so the whole expression is well behaved. 
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Series expansions of Li n (z), n > 2 . The derivations of the series expan- 
sions of Li 2 (z) presented in the previous section are by no means restricted 
to the weight two case, but we can repeat exactly the same steps for classical 
polylogarithms of arbitrary weight. In Ref. [22] the following more general 
version of Eq. (I4.2ip was proven, 

i oo An) 

u ^ z ) = "r^w logn ^ z log (~ log z ^ + £ tt logfc z ■ ( 4 - 27 ) 



The proof goes by recursion in the weight, and we refer to Ref. [35] for more 
details. 

The generalization of Eq. (I4.20p to arbitrary weight is simply given by, 

OO 1 

U n(z) = £ [B * (B®__®B®B)] k (- log(l - z)) k+l . (4.28) 

fc=0 ra-2 times 

The proof goes by recursion in the weight. We have shown in the previous 
section that Eq. (I4.28P is true for n = 2. If we assume in addition Eq. ( 14 . 2 8 [) 
true up to weight n, then we obtain, 

i- 6 "* dt 



LWi(l-0 = / -Li n (t) 
Jo t 



' _u n-2 times 



oo 



1 oo R „Z+m+l (4.29) 

n y w ji ^ ^ l + m+1 

( - u n-1 times m ~ u 

°° 1 r X k+1 



£ - [B * (S®___®_B_®S)] re 



fc— n— 1 times 



The series expansions in log z and log(l — z) for the classical polylogarithms 
are enough to evaluate all basis functions, except B^ \ for i G {16,17,18}, 
numerically in a fast and reliable way. For B^\ for i G {16, 17, 18} we need to 
extend these series expansions beyond the case of classical polylogarithms. 



4-3 Series expansions of H '(0,1,0, — 1; z) and H(0, 1, 1, — 1; z) 



In this section we present the analogues of the series expansions in log(l — z) 
and log z of the previous section for the three remaining basis functions that 
cannot be expressed in terms of classical polylogarithms only, namely B^ , for 
% G {16, 17, 18}. As the derivation of the expansions follows exactly the same 
lines as for the classical polylogarithms, we content ourselves to present the 
results and refer to Appendix[B]for details. 
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Let us start by deriving the expansion of B^\ i G {16, 17, 18} into a series in 
log(l — z). We can obtain these expansions as a corollary of a more general 
result: If a function %{z) admits a Taylor expansion of the form 

WU-^Hg^^S (4.30) 
then for a G { — 1,0,1} the functions H a {z) defined by 

n a (z) = f dtf(a;z)H(t), (4.31) 
Jo 

where f(a; z) was defined in Eq. ( 12. 2p . admit the following Taylor expansions 



n+2 



n= (^ + l)! n= (^ + 2)! 

n (l-e-n=± ( -^^^\ (4.32) 



w-i(i-c-) = £ 



(€*s(/i))r 



nt^o (n + 1)! 



■ X 



n+1 



To prove these identities we start from the integral representation (I4.3ip and, 
after having performed the change of variable t = 1 — e~* , we insert the 
series expansion for the integrand and integrate term by term, which yields 
immediately Eq. ( 14. 32ft . Combined with the expansions for HPL's of weight 
one, Eq. (14.321) allows us in principle to recursively expand HPL's of arbi- 
trary weight into a series in log(l — z). A similar iterative procedure was 
already described in Ref. [33]. As the classical polylogarithms are just a spe- 
cial case of the harmonic polylogarithms, Eq. (14. 28ft can be seen as a spe- 
cial solution of Eq. (14.321) for if(0 n _i, 1; 1 — e~ x ). Moreover, we can easily 
read off from Eq. (14.321) the corresponding expansions of H (0, 1, 0, — 1; z) and 
ff(0, 1,1,-1; 2), 



H{Q, 1, 0, -1; z) = t [B *^*f )]n (- log(l - z)Y +1 



00 [B*s 2 



(4.33) 



H(0, 1, 1, -1; z) = E 1 /I V, JW (- log(l - ^))" +1 



with s 2 (£) n = £ n -2/(n + l) . The expansions (14.331) are used inside Chaplin to 
evaluate the basis functions B^\z), i G {16, 17, 18} for z close to the origin. 
Similar to the case of the classical polylogarithms, these series converge rather 
slowly if z is close to 1. We therefore need additional expansions in log z that 
have a good convergence behavior in that region. 

In an annulus inside the unit disc, the basis functions B^\z), i G {16, 17, 18} 
can be expanded into a series in logz. As these basis functions are entirely 
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expressed through H(0, 1, 0, — 1; z) and H(0, 1, 1, — 1; z), it is enough to find 
the expansions for these cases. We obtain, 



(4.34) 



(1 \ K 1 7 4 1 
- ] - -Cs log2 + — log 4 2 
2/ 2 s3 & 480 6 & 

+ log 2 2 - - Cs log* + — log2 log 2; + — Li 2 (z) 
080 12 

+ log 2 log z U 2 (z) - 2 log 2 Li 3 (*) + - £ ^ log n+2 * , 

Z n=0 l n "I" ^J" 

#(0,1,1,-1;*) = Li 4 f-^ - - Cs log 2 + — + — log 4 2 
v ' ' ' ' ; \2J 8 s & 720 24 & 

+ ^ log 2 2 - "Q\ U 2 {z) + {I C 3 + i log 3 2 - ^ log 2) log z 
2 £0 (^ + !) ! n 



with 

7 n = — and /3 n = . (4.35) 

n n 

The proof of the Eq. (I4.34p is sketched in Appendix [Bj A similar expansion in 
the particular case of H(l, n , 1; z) has already been considered in Ref. [37] . 
We stress that the relations ( 14. 34ft are only valid in the region Re(z) > and 
R < \z\ < 1, for some R (we give a precise value for R in the next section 
when discussing the implementation of these expansions into Chaplin). For 
Re(z) < 0, we proceed in following way, 

• for H(0, 1,0, —1; z), we map the problem to Re(z) > via the functional 
equation 

#(0,1,0,-1;*) = -#(0, -1,0,1;-*) 

1 - ■ /■ 2 x 1- ( z 2 



#(0, 1, 0, -1; -z) - - Li 4 1 - z 2 ) + -Li 4 



2 ^ j 2 V* 2 - 1 

- 2 Li 4 (^— \ + 2 UA-z) + 2 LiAz) - 2 Li 4 , 

VI — zj \z + 1, 

- 2 Li 4 ( j^y) + 2 Li 4 (l + z)-2 U 3 (-z) log(l - z) (4.36) 

- 2 Li 3 (*) log(l + z) - I Cs log(l - z) - I Cs log(l + z) 

- ^ log 4 (l - *) - 1 log 4 (l + 2) + i log(-z) log 3 (l - z) 
+ — log(l + z) log 3 (l - *) + — log 3 (l + z) log(l - z) 
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+ ^ log(-z) log 3 (l + 2) + ^ log 2 (l + ^) log 2 (l - z) 

1 57T 2 

- - log(-z) log(l + z) log 2 (l -z) + — log 2 (l - z) 

I 24 (43?) 

- - log(-z) log 2 (l + z) log(l - *) - — log 2 (l + z) 

2 4 

- U 2 {-z) Li 2 (z) + ^ log(l + z) log(l - z ) + ^. 

• for iJ(0, 1, 1, — 1; z), we use the following expansion, valid for Re (2) < 0, 

H(0, 1, 1. -1; *) = 2 Li 4 - — + — log 4 2 + — log 2 2 
+ Q Cs - I log 3 2 + log2) log 2:- log 2 2^ Li 2 (-z) 

- 4 5 7^TT)T n log { ~ z) l log( - log( ^ )} - ^TT~nJ 

1 - (l7*,( 2 )) w 1 +1 1 g ^*(g«Q] w j w+2( _, 

+ 4 n ^ (n + 1)! n° g 1 " J + 4^ („ + 2)! l0g 1 " j ' 

(4.38) 

with 

g W = ^ an d £ n = C-n • (4.39) 

At this stage we have all the ingredients to evaluate all harmonic polyloga- 
rithms up to weight four numerically for arbitrary complex arguments. To this 
effect, we have implemented the decomposition of HPL's to the basis defined 
by Eq. (13.31 - 13.9ft and the series expansions presented in this section into the 
Chaplin library, which we will present in the next section. 



5 The Fortran library CHAPLIN 



5. 1 Installation 



The Chaplin code is available as a .tar archive on the website 
htt p : //projects .hepforge . org/ chaplin/| . 

Once unpacked, the code can be compiled via 

. /configure 
make install 
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As a result, both a static and a shared library are created in the directory 
/usr/lib, which may require root privileges. The directory which the library 
is installed in can be changed during configuration via 

./configure — pref ix=/"path to Chaplin/" 

This feature allows to create the library even without root privileges. 

We advocate the use of the shared library when linking CHAPLIN to a Fortran 
code, such that the routines needed by the program will only be called dur- 
ing run-time. Static linking, on the contrary, puts all the Chaplin functions 
into the executable Chaplin is linked against, which might result in long 
compilation times and rather large executables. 

5.2 Running CHAPLIN 

Once compiled, the Chaplin library can be linked to other programs in the 
same way as any other library. This enables the user to call the numerical 
routines of CHAPLIN from within his/her own program. The function calls to 
the numerical routines of Chaplin are done via the following functions, 

double complex HPLl(nl, z) 

double complex HPL2(nl, n2, z) 

double complex HPL3(nl, n2, n3, z) 

double complex HPL4(nl, n2, n3, n4, z) 

with ni G {—1,0,1} and z being any double complex number. The return 
value is H(n; z), i.e., the HPL with weight vector n at the point z. Alterna- 
tively, in case only the real and/or imaginary parts of an HPL are needed, the 
user might find the following functions convenient, 

double precision HPL2real(nl, n2, xr, xi) 
double precision HPL2im(nl, n2, xr, xi) 

and similarly for HPL3, HPL4. The variables xr, xi are double precision 
variables, denoting the real and imaginary parts of the argument of the HPL. 
Note that these functions are useful to call Chaplin from within a C++ pro- 
gram, where complex numbers are not natively supported. 

When one of the aforementioned functions is called, Chaplin starts by de- 
composing the corresponding HPL internally into the basis of 32 functions 
described in Section [3j In a second step, Chaplin proceeds to the numerical 
evaluation of the individual basis functions appearing in the decomposition. 
The numerical routines called to this effect are different depending on the 
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Fig. 1. The different regions of the complex plane used inside the Chaplin library. 

value of the argument z of the HPL. More precisely, the complex plane is 
divided into six regions (See Fig. [TJ, 

• Region I: inside an annulus 0.025 < \z\ < 0.3, the basis functions are eval- 
uated by using the expansions in log(l — z) presented in Section HI 

• Region II: inside an annulus 0.3 < \z\ < 1, the basis functions are evaluated 
by using the expansions in log z presented in Section HJ 

• Region III: points outside the unit disc, \z\ > 1, are mapped back to the 
interior of the unit disc via inversion relations. 

• Regions IV & V: The basis of Section [3] involves functions that are loga- 
rithmically divergent at ±1, leading to spurious singularities in the basis 
expansion. To avoid numerical instabilities caused by these spurious singu- 
larities, we use Taylor expansions close to z — ±1 to evaluate the individual 
HPL's without proceeding to a decomposition into basis functions. 

• Region VI: In order to achieve a good numerical precision close to the origin 
of the complex plane, we use Taylor expansions in a disc \z\ < 0.025 without 
proceeding to a decomposition into basis functions. 

At the end of this procedure, Chaplin returns the numerical value of the HPL 
given as an input. In case the numerical evaluation of a divergent quantity 
is attempted (e.g., H(0 n ; 0) or H(±l, a; ±1)) an exception is thrown and the 
evaluation is aborted. Note that for real values of the argument, Chaplin uses 
the prescription, i.e., for zel, H(a; z) is interpreted as H(a; z + ie). 

We conclude this section by giving an example of a sample program that prints 
all nine HPL's of weight two at the point z = 1.54 + 0.91i. 
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program chaplintest 

double complex HPL2,z 
integer nl,n2 

z = dcmplx(1.54d0,0.91d0) 
do nl=-l,l 

do n2=-l,l 

print*, HPL2(nl,n2,z) 

enddo 
enddo 

end program 

5.3 Validation and comparison to other codes 

We have compared the results obtained by Chaplin against hplog and GiNaC 
for real values of the argument and against HPL and GiNaC for complex 
values of the argument. The results for some sample points are summarized 
in Tab. [2]- |5j For real arguments, a small imaginary part is assumed. Note 
that the reduction of H{1, —1, — 1, 0; z) to basis functions is among the most 
complicated cases, involving almost 50 lines of Fortran code. 

In order to give an idea of the CPU time needed per function call, we present 
in Tab. [1] the time for one million function calls on a 3 GHz Intel X5450 
Processor for some HPL's in the six different regimes of the complex plane 
shown in Fig. [TJ The running time varies only marginally inside a given region. 
The resulting average times for a single function call are given in Tab. [TJ in 
units of microseconds {us). 



Region 


I 


II 


III 


IV 


V 


IV 


Li 2 (z) 


1.4 


4.4 


4.5 


0.2 


0.2 


0.3 


Li 3 (z) 


1.9 


4.7 


5.0 


0.2 


0.2 


0.3 


hu{z) 


6.2 


8.9 


9.1 


4.2 


4.3 


4.5 


#(0,1, 0,-1; 2:) 


8.7 


17.3 


44.2 


4.3 


4.4 


4.5 


ff(l,-l,-l,0;z) 


140.0 


213.4 


393.8 


4.3 


4.8 


4.7 



Table 1 

Average time per function call in microseconds (/us), for the six regions of the 
complex plane shown in Fig. [TJ 
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z = 0.5 




Chaplin 
hplog 
GiNaC 


5.8224052646501201e-01 + 0. i 
5.8224052646501256e-01 + 0. i 
5.8224052646501245e-01 + 0. % 




Chaplin 
hplog 
GiNaC 


5.3721319360804021e-01 + 0. % 
5.37213193o0804010e-Ul + 0. i 
5.3721319360804021e-01 + 0. i 




Chaplin 
hplog 
GiNaC 


5.1747906167389901e-01 + 0. % 
5.174790olb7389945e-01 + 0. i 
5.1747906167389934e-01 + 0. i 


H(0, 1,0, -l;z) 


Chaplin 
hplog 
GiNaC 


7.7856141848313908e-02 + 0. i 
/ . ^ 85bl4184o31o090e-02 + 0. z 
7.7856141848313215e-02 + 0. i 


H(l, -1,-1,0;*) 


Chaplin 
hplog 
GiNaC 


-6.3908284909225732e-02 + 0. i 
-6.3908284909226009e-02 + 0. i 
-6.3908284909226135e-02 + 0. i 



Table 2 

Comparison between Chaplin, hplog and GiNaC for the point z = 0.5. 
6 Summary 

In this paper we have presented CHAPLIN, a new Fortran library to compute 
harmonic polylogarithms up to weight four for arbitrary complex argument. 
The algorithm is based on a reduction of HPL's to a set of basis functions 
which are then evaluated numerically using series expansions allowing for a 
very fast numerical convergence, hence rendering the computational cost of 
a function call quite modest. We have checked our numerical results against 
well-established codes [3D]I31|I32|[H5] and found agreement to at least 14 digits 
for any argument in the complex plane. 
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z = 2.0 


Li 2 (z) 


Chaplin 
hplog 

GiNaC 


2.46740ll002723404e+00 + 2.l775860903036022e+00 i 
2.46740ll002723399e+00 + 2.l7758609030360l7e+00 i 
2.46740ll002723395e+00 - 2.l775860903036022e+00 i 


Li 3 (^) 


Chaplin 
hplog 

GiNaC 


2.76207l9062289245e+00 + 7.5469382946024677e-0l i 
2.76207l906228924le+00 + 7.5469382946024799e-0l i 
2.76207l906228924le+00 - 7.54693829460248l0e-0l i 


Li 4 (z) 


Chaplin 
hplog 

GiNaC 


2.4278628067547032e+00 + l.7437l30002545320e-0l i 
2.4278628067547032e+00 + 1.7437130002545298e-01 % 
2.4278628067547032e+00 - 1.7437130002545306e-01 % 


H(0,l,0,-l;s) 


Chaplin 
hplog 

GiNaC 


5.1994752047739468e-01 + 1.7909927717176168e+00 i 
5.1994752047739512e-01 + 1.7909927717176164e+00 i 
5.1994752047739445e-01 - 1.7909927717176168e+00 i 


H(l,-l,-l,0;z) 


Chaplin 
hplog 

GiNaC 


8.0548200591357266e-01 - 1.3189461296972327e+00 i 
8.0548200591356789e-01 - 1.3189461296972333e+00 i 
8.0548200591356811e-01 + 1.3189461296972318e+00 i 



Table 3 



Comparison between Chaplin, hplog and GiNaC for the point z = 2.0. Note that 
GiNaC uses a different convention for the imaginary parts. 
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z = 0.5 + 0.5 i 




Chaplin 
HPL 
GiNaC 


4.53985269l5029508e-0l + 6.4376733288926902e-0l i 
4.53985269l502954le-0l + 6.4376733288926879e-0l i 
4.53985269l5029558e-0l + 6.4376733288926880e-0l i 


Lis(«) 


Chaplin 
HPL 
GiNaC 


4.86l5953708555987e-0l + 5.7007740708876864e-0l i 
4.86l59537085560l4e-0l + 5.7007740708876930e-0l i 
4.86l5953708556009e-0l + 5.7007740708876897e-0l i 


Li 4 (z) 


Chaplin 
HPL 
GiNaC 


4.9578ll2l82l8388le-0l + 5.3402238407975344e-0l i 
4.9578ll2l82l83897e-0l + 5.3402238407975377e-0l i 
4.9578ll2l82l83876e-0l + 5.3402238407975355e-0l i 


H(0, l, 0, — l; z) 


Chaplin 
HPL 
GiNaC 


-3.6325772l79994l09e-02 + l.38499l682646747e-0l i 
-3.6325772179994879e-02 + 1.3849916826467456e-01 i 
-3.6325772179994845e-02 + 1.3849916826467457e-01 i 


H(l,-l,-l,0;*) 


Chaplin 
HPL 
GiNaC 


9.1142643382278842e-02 - 9.8191320890700317e-02 i 
9.1142643382278176e-02 - 9.8191320890700678e-02 i 
9.1142643382278163e-02 - 9.8191320890700595e-02 i 



Table 4 

Comparison between Chaplin, HPL and GiNaC for the point z = 0.5 + 0.51 
A Proof of the associativity of the convolution product 



In this appendix we proof the associativity of the convolution product, Eq. (14. 2p . 
Using the definition of the convolution product, Eq. (14. ip . we obtain, 

((a * b) * c) N = Y (] ( a * b )m c N - m = Y Y ( J [j a n b m - n c N „ m . 

m=0 \ m J m=0n=0 \ m J \ H ) 

(A.l) 

We now exchange the sums over m and n and shift the summation variable 
such that all sums start from zero. This gives, 



Q"n bm—n C/V- 



» N -( N \(m + n\ 

VV \a n b m c N _ m _ n (A.2) 

^ Q ^ \ m + n J\ n J 

" » / N \fm + N-n\ 
kL\m + N-n)[ N-n )^h m ^ m , 
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z = 2.0 + 2. 0i 




Chaplin 
HPL 
GiNaC 


3.44973l2626l78323e-0l + 2.7342872l86403562e+00 i 
3.44973l2626l78264e-0l + 2.7342872l86403560e+00 i 
3.44973l2626l7826le-0l + 2.7342872l86403562e+00 i 


Lis(«) 


Chaplin 
HPL 
GiNaC 


l.237054890750l702e+00 + 2.70246078223l0069e+00 i 
l.237054890750l697e+00 + 2.70246078223l0064e+00 i 
l.237054890750l697e+00 + 2.70246078223l0065e+00 i 


Li 4 {z) 


Chaplin 
HPL 
GiNaC 


l.7008027579027265e+00 + 2.4625762l77390942e+00 i 
l.7008027579027260e+00 + 2.4625762l77390939e+00 i 
l.700802757902726le+00 + 2.4625762l77390937e+00 i 


H(0, l, 0, -l; z) 


Chaplin 
HPL 
GiNaC 


-l.309292l033357463e+00 + 8.60095l353690l472e-0l i 
-1.3092921033357458e+00 + 8.6009513536901561e-01 % 
-1.3092921033357459e+00 + 8.6009513536901561e-01 i 


H(l,-l,-l,0;*) 


Chaplin 
HPL 
GiNaC 


1.3154184588794104e+00 - 2.6274818437873471e-01 % 
1.3154184588794054e+00 - 2.6274818437872689e-01 % 
1.3154184588794056e+00 - 2.6274818437872688e-01 i 



Table 5 

Comparison between Chaplin, HPL and GiNaC for the point z = 2.0 + 2.0 i. 

where the last step follows from changing the summation variable according 
to n — > N — n. The product of binomials can be simplified, 



' N 
,m + N 



71, 



f m + iV — n 
. N-n , 



(N + m- n)\ 



(N + m — n)\(n — m)\ m\(N 
N\ n\ (N\ 



n\{N — n) m\(n — m)\ \n 



■ n) 
m, 



(A.3) 



yielding, 

((a * b) * c) N = E 



N n 



n=0 m=0 



N 

, n , 



n 



(a*(b*c)) N . (A.4) 



B Derivation of Eq. (Olj) 



In this appendix we present the derivation of the series expansions given in 
Eq. (I4.34p . The derivations follow the same spirit as in all other cases discussed 
in Section HI i.e., we start from the integral representation and perform a 
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change of variable, before inserting the Taylor expansions of the integrand. 
The proof involves some technical issues, which are discussed in detail in the 
following on the example of H(0, 1, 0, — 1; e x ). All other cases are similar. 

To derive the Taylor expansion of H(0, 1, 0, — 1; e x ), we start from the series ex- 
pansion of H(—l; e x ) = log(l + e x ) and then integrate up to H(0, 1, 0, — 1; e x ). 
The series expansion of H(—l;e x ) = log(l + e x ) is easily obtained from the 
integral representation and the generating function of the Genocchi numbers, 
Eq. (I4.13p . Letting t = e*', we get 



i_l + t & 2 Jo f 1 + e_-*' 

oo 



1 ™ 'V 

log2-- V ^a;™ 
z n=o ' i - 



where we introduced the shorthand 



7n = — • (B.2) 



n 



Next we determine the series expansion of H(0, — lj eX ) — — Li 2 (— e x ). Repeat- 
ing exactly the same steps as for H(—l; e x ), we obtain 



re x Af :~ - fX 

H(0, -1; e") = -Li a (-1) + / ^ H(-l; t) = - + / dt' (-1; e*) 

ji t 12 Jo 



12 2^(n + l)! 



(B-3) 



If we try to repeat the same procedure for H(l, 0, — 1; e x ), a technical difficulty 
arises: in the previous example we had the split the integral over [0, e x ] into two 
contributions from [0,1] and [l^]. In the present case, however, we cannot 
do this, because H(l, 0, — 1; 1) is divergent. We therefore introduce a regulator 
e and split the integration region into [0, e £ ] and [e £ , e x ], and we will take the 
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limit £ — > at the very end. We then obtain, 

r eX dt 

ff (1, 0, -1; e x ) = ff (1, 0, -1; e £ ) + / -— ff(0, -1; t) 

= ff (1, 0, -1; e^) - f H(0, -1; e*') 

= ^(l,0,-l;e B ) + ^[^(l;e')-^(l;^)] 



+ log 2 [if (1, 0; e x ) - ff (1, 0; e £ )] + - £ r 



1 « (g * 7)n , 

(B.4) 



We now have to send the regulator to zero. Concentrating only on terms 
depending on e and using the shuffle algebra for HPL's, we obtain, 

ff (1, 0, -1; e £ ) - ^ ff (1; e £ ) - log 2 #(1, 0; e e ) 

= ff(l;e £ )ff(0,-l;e £ )-^ff(l;e e ) (R5) 
- if (0, 1, -1; e £ ) - ff (0, -1, 1; e £ ) - log 2 ff (1, 0; e £ ) . 

All the terms in the second line have a smooth limit as e — > 0, and the two 
divergent terms in the first line cancel exactly. This leaves us with 

ff (1, 0, -1; e x ) = ~ Ca - ^ log 2 + ^ ff (1; e») + log 2 ff (1, 0; e*) 

1 ~ (fj* 7 ) w r „ +1 (R6) 



The last integration is easy to perform and we immediately obtain 

(1 \ T 7 4 "i 
- J - -Ca log2 + — log 4 2 
2/ 2 s3 & 480 6 & 

2 r 2 2 

+ 1- log 2 2 - - C 3 x + y log 2 x + ^ Li 2 (e*) (B.7) 



+ log 2 x Li 2 (e x ) - 2 log 2 Li 3 (e 3; ) + - £ 



1 ^ (S*7) ? 



2 n =o(^ + 2)! 

which agrees with Eq. fl4.34j) after replacing x = logz. 



x 



n+2 
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